{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to generate rates of elevation change and time series from ICESat-2\n",
    "\n",
    "ICESat-2 Hackweek 2020\\\n",
    "Johan Nilsson\\\n",
    "Jet Propulsion Laboratory\\\n",
    "2020-06-16\n",
    "\n",
    "In this tutorial we will present an easy processing flow to generate gridded altimetry time series of elevation change based on the captoolkit software. We will then use these results to produce estimates mass balance for Pine Island Glacier. Its should be noted that this processing flow works equally well for other altimetry missions as well, both space- and airborne.\n",
    "\n",
    "## Tutorial preparation\n",
    "\n",
    "1. Log on to the event Jupyterhub: https://icesat-2.hackweek.io\n",
    "2. Open a terminal in Jupyterlab\n",
    "3. Clone repo: git clone https://github.com/ICESAT-2HackWeek/ScienceDataGeneration.git\n",
    "4. Enter the directory and open the notebook\n",
    "\n",
    "## Introduction\n",
    "\n",
    "In many aspects of Earth Science gridded datasets plays an important role for understanding and extracting geophysical signals in both space and time. \n",
    "\n",
    "So what are the reasons why we would like to work with gridded datasets?\n",
    "\n",
    "- Understaing geophysical signals in both time and space based on noisy observations \n",
    "- Reducing data to a size that is usable for storage \n",
    "- To allow for easy analysis and visualization of large scale data sets  \n",
    "- Application of advanced digital signal processing algorithms.\n",
    "\n",
    "Satellite altimetry has for decades been used to study the worlds ice sheets and glaciers by measuring the change of elevation over time. We will here present a method to separate static topography from the time variable signal of interest. This will allow us to both create gridded datasets for time-variable elevation change and to generate DEMs bases on the static topography. \n",
    "\n",
    "## Learning objectives\n",
    "- Get familiar with a simple altimetry processing flow to generate elevation change and time series\n",
    "- Understand what the important parameters are in the flow\n",
    "- Understand limitations and bottleneckz\n",
    "- Tuneing of parameters obtain good results \n",
    "- Gridding and filtering off different geophysical entities\n",
    "\n",
    "![is_dhdt.jpeg](is_dhdt.jpeg \"Figure of Eleavtion Change from ICESat\")\n",
    "[ Image credit: Johan Nilsson (JPL) ]\n",
    "\n",
    "## Data products and general approach\n",
    "\n",
    "The data we will use in this excercis is the ATL06 data product, which has been flatten and sub-sampled based on the previous excercies. It contains the minumim number of variabels needed to run this example and generate sciece worthy data products. In this case we only need: **lon, lat, time and height**.\n",
    "\n",
    "We will here show you how to make the following a simple way to make a few intresting and important science products:\n",
    "\n",
    "- Time Sereis of Elevation Change\n",
    "- Computing Rates of Elevation Change\n",
    "- Constructing high-resolution DEMs\n",
    "- Quality control of altimetry data\n",
    "\n",
    "These are already standard products from the ICESat-2 mission which means you don't have to do it yourself...nice right! Big thanks to Ben Smith for putting all the hard work to make that happen! However, the processing flow presented here is generic and can be applied to any altimeter mission."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Case study: Pine Island Glacier - Cookbook recipe\n",
    "In this section we will use surface elevations from ICESat-2 from 2018/10 - 2020/04 over Pine Island Glacier (PIG) which is one the most active glaciers in Antartcica. But before we start with the excersice we need to import some packages and functions that we need for the analysis. Many of them can be found in the \"utils.py\" file and have been written and devloped here at JPL as part of the **captoolkit**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import numpy as np\n",
    "import cartopy.crs as ccrs\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from utils import transform_coord\n",
    "from utils import make_grid\n",
    "from utils import mad_std\n",
    "from utils import spatial_filter\n",
    "from utils import interp2d\n",
    "from utils import tiffread\n",
    "from utils import binning\n",
    "from scipy.ndimage.filters import generic_filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, before we do anything lets have a look at the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ../../tutorial-data/gridding-time-series/org/IS2*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You notice that we have two different files with the endings \"A\" and \"D\". This is due to the fact that we have separated the data into ascending and descending orbits. This was done directly using captoolkit's ATL06 reader **readatl06.py**, which is just one of many readers avalibale. The need for this will become clear very soon I hope when we start processing the data.\n",
    "\n",
    "Let's look at the variables in the data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!h5ls ../../tutorial-data/gridding-time-series/org/IS2_HEIGHTS_PIG_A.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now when we have done that let's load the ICESat-2 ATL06 elevations into memory and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('../../tutorial-data/gridding-time-series/org/IS2_HEIGHTS_PIG_A.h5','r') as f_a:\n",
    "    lat   = f_a['lat'][:]\n",
    "    lon   = f_a['lon'][:]\n",
    "    h_elv = f_a['h_elv'][:]\n",
    "    t_yrs = f_a['t_year'][:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also read in a grounded ice mask that we can use later to edit our data with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('../../tutorial-data/gridding-time-series/ANT_groundedice_1000m.h5','r') as f_m:\n",
    "    Xm = f_m['X'][:]\n",
    "    Ym = f_m['Y'][:]\n",
    "    Zm = f_m['Z'][:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the data and see what we get!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(9,9))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "ax.set_extent([-180, 180, -90, -66], crs=ccrs.PlateCarree())\n",
    "plt.scatter(lon[::100], lat[::100], s=3, c=h_elv[::100], alpha=.7, transform=ccrs.PlateCarree(), cmap='terrain')\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation (m)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can here see our region of interest (ROI) that we will use for the analysis. Now will zoom in so we can get more detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(9,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.scatter(lon[::100], lat[::100], s=3, c=h_elv[::100], alpha=.7, transform=ccrs.PlateCarree(), cmap='terrain')\n",
    "plt.clim([100,1400])\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation (m)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, you can see how data density of ICESat-2 tracks looks over Pine Island. You can for example see we have some data missing here and there due to intial filtering in the reading procedure (quality summary variable mostly). We can clearly see trunk and the catchment of PIG here just looking at the elevation. Remeber we are only plotting every 100 points so there is 100 times more data than shown here!!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Min. height:', int(h_elv.min()),'(m)', 'Max. height:', int(h_elv.max()),'(m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cm = plt.cm.get_cmap('terrain')\n",
    "Y,X = np.histogram(h_elv, 70)\n",
    "x_span = h_elv.max()-h_elv.min()\n",
    "C = [cm(((x-X.min())/x_span)) for x in X]\n",
    "plt.bar(X[:-1],Y,color=C,width=X[1]-X[0],edgecolor='k')\n",
    "plt.ylabel('Frequency')\n",
    "plt.xlabel('Elevation (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## captoolkit - Cryosphere Altimetry Processing Toolkit\n",
    "\n",
    "https://github.com/fspaolo/captoolkit\n",
    "\n",
    "This repo is a collection of command-line programs to read, process and analyze altimetry datasets. It is a consolidation of several years of work and experience in how to work with altimetry data to generate science data to understand our changing Cryosphere. \n",
    "\n",
    "Here, we are going to show how to use a small portion of these programs to create a simple processing workflow that will provide us the ability to go from elevation point data point to gridded data products of time series and elevation. This software is used daily for our own research and is constantly changing and hopefully improving, so please let us know if we can make any improvement or if something is not working properly.\n",
    "\n",
    "In this excercies we are going to use four different programs avalibale in the repo:\n",
    "\n",
    "- **fittopo.py**: Detrend surface elevation to separate time-variable elevation from static topography by removing a matehmatical surface.\n",
    "- **merge.py**: simple program for merging hdf5-files and selecting variables wanted variables.\n",
    "- **interpgaus3d.py**: Spatio-Temporal interpoaltion of irregular data to grid using gaussian weights.\n",
    "- **readatl06.py**: reading and separation of ATL06 into induvidual beams files\n",
    "- **xover.py**: program for computing satellite crossover differences \n",
    "\n",
    "That’s it! Those are the main program that you need to construct time series! But we have also gathered plenty of helpful functions in **utils.py** that can be imported for further data filter, gridding and analysis. We will show some of these here in the notebook and how you can use them.\n",
    "\n",
    "Here is a summary for the functions used from **utils.py**:\n",
    "\n",
    "- **make_grid**: construct a grid provided boundries and resolution\n",
    "- **binning**: binning of data using mean or median with corresponding statistics \n",
    "- **spatial_filter**: spatial outlier filter using robust n-sigma rule and threshold \n",
    "- **transform_coord**: coordinate transformation using pyproj based on EPSG-number\n",
    "- **interp2d**: point to raster interpolation based on ndimage's map_coordinates\n",
    "- **interpkrig**: interpolation using ordinary kriging/collocation of irregular data in space\n",
    "- **interpmed**: interpolation using median of irregular data in space\n",
    "- **interpgaus**: interpolation using gaussian kernel of irregular data in space\n",
    "- **tiffread**: reading of tiff-files\n",
    "\n",
    "Here is the link to the repo if you want to download and use it later on your own personal computer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#git clone https://github.com/fspaolo/captoolkit.git\n",
    "#cd captoolkit\n",
    "#python setup.py install"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok ones you have installed the repo let's have a look at the first program **fittopo.py** which gives you the opportunity to separate the time-variable surface elevation form the static one:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python fittopo.py -h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Simplified this program overlays a grid of desired resolution on top of the data and fits to the data around each grid-node, given a provided search radius, either with a biquadratic (1), bilinear (2) or \"mean\" (3) model. To improve the quality of the solution the program uses a random-walk approach to find the centroid of the data inside the search radius to \"lock\" on to repeat-tracks for example.\n",
    "\n",
    "We are going to run ascending and descending tracks separately this is good practice as mixing them both in one solution can introduce errors as there might be A/D bias present. In our topographical model we also include a linear rate to make sure that all data are centered to a common time. This has another positive effect as it allows us to use all the mean heights from the solution to generate a DEM, as they are now centered to the same time.\n",
    "\n",
    "Let us have a look at program and the input!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./fittopo.py ../../tutorial-data/gridding-time-series/org/IS2_HEIGHTS_PIG_A.h5 -d 1 1 -r 1.0 -q 3 -i 5 -z 10 -m 100 \\\n",
    "-k 2 -t 2020 -j 3031 -v lon lat t_year h_elv -s 10 -p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!fittopo.py ../../tutorial-data/gridding-time-series/org/IS2_HEIGHTS_PIG_D.h5 -d 1 1 -r 1.0 -q 3 -i 5 -z 10 -m 100 \\\n",
    "-k 2 -t 2020 -j 3031 -v lon lat t_year h_elv -s 10 -p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code using 1 km grid spaceing takes about ~7 min to run so we have already prepped the data for you!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ../../tutorial-data/gridding-time-series/proc/IS2*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!h5ls ../../tutorial-data/gridding-time-series/proc/IS2_HEIGHTS_PIG_A_TOPO.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should note now that we have added an extra set of variables to the file! For us **h_res** is the most interesting as it contains the time variable residuals. We also have **e_res** which is the RMSE from the fit and can be used to weight our data. **h_mod** is the intercept or mean height from the solution and is very useful to create DEMs and **slp** contain the surface slopes in the x- and y direction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merging of files\n",
    "Now we have removed the static topography from our ascending and decending orbits and thus we need to merge the two files. captoolkit come with such capabilities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./merge.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./merge.py ../../tutorial-data/gridding-time-series/proc/IS2_HEIGHTS_PIG_*.h5 -o ./IS2_HEIGHTS_PIG_AD.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we created a new file in our the directory by merging the two files, lets see! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls *.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del(lon,lat) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's read this file and try to construct some time series!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('./IS2_HEIGHTS_PIG_AD.h5','r') as f_r:\n",
    "    lat   = f_r['lat'][:]\n",
    "    lon   = f_r['lon'][:]\n",
    "    t_yrs = f_r['t_year'][:]\n",
    "    h_res = f_r['h_res'][:]\n",
    "    h_mod = f_r['h_mod'][:]\n",
    "    h_org = f_r['h_elv'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.scatter(lon[::100], lat[::100], s=3, c=h_res[::100], alpha=.7, transform=ccrs.PlateCarree(), cmap='coolwarm_r')\n",
    "plt.clim([-5,5])\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation residual (m)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.hist(h_res[~np.isnan(h_res)], 100, edgecolor='k')\n",
    "plt.xlabel('Residuals (m)')\n",
    "print('Mean:', np.around(100*np.nanmean(h_res),3), 'cm')\n",
    "print('Std.dev:',np.around(100*np.nanstd(h_res),3), 'cm')\n",
    "print('Rejected:', np.around(100*(float(len(h_res[np.isnan(h_res)]))/len(h_res)),2), '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well this looks noisy! Lets try to filter the data a bit using a global threshold of 10 m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_res[np.abs(h_res) > 10] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.hist(h_res[~np.isnan(h_res)], 100, edgecolor='k')\n",
    "plt.xlabel('Residuals (m)')\n",
    "print('Mean:', np.around(100*np.nanmean(h_res),3), 'cm')\n",
    "print('Std.dev:',np.around(100*np.nanstd(h_res),3), 'cm')\n",
    "print('Rejected:', np.around(100*(float(len(h_res[np.isnan(h_res)]))/len(h_res)),2), '%')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is a lot better now we can clearly see a normal distribution around zero. But we can maybe do a bit better! Let's apply some outlier filtering that works in the spatial domain. For this we need to functions that are available in **utils.py**. The first one is **transform_coord** which allows us to move between different projection and is wrapper for pyproj."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_i, y_i = transform_coord('4326', '3031', lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another fucntion in **utils** we now use is called **spatial_filter** which applies a robust n-sigma filter on the residuals inside a given domain. In this case 100 x 100 km with n=3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_i = spatial_filter(x_i, y_i, h_res.copy(), dx=100e3, dy=100e3, n_sigma=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Mean:', np.around(100*np.nanmean(z_i),3), 'cm/yr')\n",
    "print('Std.dev:',np.around(100*np.nanstd(z_i),3), 'cm/yr')\n",
    "print('Rejected:', np.around(100*(float(len(z_i[np.isnan(z_i)]))/len(z_i)),2), '%')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.hist(z_i[~np.isnan(z_i)],100,edgecolor='k')\n",
    "plt.xlabel('Residuals (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That clearly helped a little bit to remove some more noise! Let us plot the residuals again to loook what changed!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.scatter(lon[::100], lat[::100], s=3, c=z_i[::100], alpha=.7, transform=ccrs.PlateCarree(), cmap='coolwarm_r')\n",
    "plt.clim([-5,5])\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation residual (m)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data now looks way more clear and we can even see some patterns, but not all the noise is gone but that's ok! Now let us investigate how the temporal variability looks like by plotting the residuals as a function of time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,4))\n",
    "xb,yb = binning(t_yrs, z_i, xmin=2018.5,xmax=2020.26, dx=0.125 / 12., window=0.25 / 12., median=True)[0:2]\n",
    "plt.plot(t_yrs[::10],z_i[::10],'.',markersize=0.1)\n",
    "plt.plot(xb, yb,'-r', linewidth=3)\n",
    "plt.ylim([-4,4])\n",
    "plt.ylabel('Elevation residual (m)')\n",
    "plt.xlabel('Time (yrs)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is a trend here! PIG has been loosing elevation for the last 1.5 years!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "p = np.polyfit(t_yrs[~np.isnan(z_i)],z_i[~np.isnan(z_i)],1)\n",
    "print('Trend:',np.around(p[0],2),'(m/yr)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatio-temporal interpolation\n",
    "Ok we now have tons of data here! Remember ICESat-2 ATL06 product has an along-track sampling of 20 m, so we need to reduce this data into workable format for further analysis. One way of doing this is to interpolate the data to a regular grid with a specific resolution. In our tutorial we are going to do this in both space and time to create monthly time series in the form of a data cube. A simple program for this purpose is **interpgaus3d.py** available in captoolkit (there is many more that work both in 2d and 3d). This is a simple 3D (or you might want to call it 4D) interpolator using a gaussian kernel in both space and time to compute a value at each grid node.\n",
    "\n",
    "**Later in the tutorial we are going to discuss different interpolation techniques in more detail**\n",
    "\n",
    "Let have a look at the program and its input!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./interpgaus3d.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./interpgaus3d.py ./IS2_HEIGHTS_PIG_AD.h5 ./data_cube.h5 -d 10 10 -t 2018.8 2020.25 0.125 \\\n",
    "-r 30 -a 10 0.125 -p 3031 -c 100 3 10 -v lon lat h_res t_year dummy -s 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let us have a look at the simple structure of the data using **h5ls**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!h5ls data_cube.h5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data are saved as arrays and vectors: **Z_pred** = gridded values, **Z_rmse** = the error for each value and **Z_nobs** = number of points in solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('./data_cube.h5','r') as f_c:\n",
    "    Xi = f_c['X'][:]\n",
    "    Yi = f_c['Y'][:]\n",
    "    ti = f_c['time'][:]\n",
    "    Zi = f_c['Z_pred'][:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolate our mask to the extent of our prediction grid using **interp2d** using the nearest point (order=0 and for bilinear order=1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zmi = interp2d(Xm, Ym, Zm, Xi.ravel(), Yi.ravel(), order=0).reshape(Xi.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lon, Lat = transform_coord('3031','4326', Xi, Yi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the mask to our datacube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zi[:, Zmi == 0] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig = plt.figure(figsize=(9,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.pcolormesh(Lon, Lat, Zi[130,::], transform=ccrs.PlateCarree(), cmap='coolwarm_r')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Quick way of constructing a ROI time series of the cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hi = np.nanmean(np.nanmean(Zi,1),1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(ti, hi,'r',linewidth=3,label='interpolated')\n",
    "plt.plot(t_yrs[::100], h_res[::100], '.', markersize=0.3,linewidth=3)\n",
    "p1 = np.polyfit(ti[~np.isnan(hi)], hi[~np.isnan(hi)],1)\n",
    "p2 = np.polyfit(t_yrs[~np.isnan(z_i)], z_i[~np.isnan(z_i)],1)\n",
    "plt.plot(ti, np.polyval(p2, ti),'k',label='linear-fit')\n",
    "plt.legend()\n",
    "print('Interpolated time series:',np.around(p1[0],2))\n",
    "print('Linear fit to data:',np.around(p2[0],2))\n",
    "plt.ylim([-3,3])\n",
    "plt.xlabel('Time (yrs)')\n",
    "plt.ylabel('Elevation residuals (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(ti, hi - hi[0],'r',linewidth=2)\n",
    "plt.ylabel('Elevation Change (m)')\n",
    "plt.xlabel('Time (yrs)')\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok we checked overall time series what about the spatial pattern of elevation change? We can get that by fitting a linear rate to each grid-cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nt,nx,ny = Zi.shape\n",
    "rate = np.zeros((nx,ny))*np.nan\n",
    "for i in range(nx):\n",
    "    for j in range(ny):\n",
    "        dh = Zi[:,i,j]\n",
    "        if np.any(np.isnan(dh)):continue\n",
    "        rate[i,j] = np.polyfit(ti,Zi[:,i,j],1)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ones we estimated the rate let's apply some gentle smoothing using a mean filter kernel. You can use any function here if you want! Also, lets use the mask that we read in before to mask out unwanted areas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zf = generic_filter(rate.copy(), np.nanmean, 3)\n",
    "Zf[Zmi == 0] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "im = plt.pcolormesh(Lon, Lat, Zf, transform=ccrs.PlateCarree(), cmap='coolwarm_r',vmin=-1.5,vmax=0.5)\n",
    "ax.coastlines('50m')\n",
    "plt.colorbar(label='Elevation Change (m/yr)')\n",
    "print('Trend:',np.around(np.nanmean(Zf),2),'(m/yr)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dVdt = np.nansum(Zf)*10e3*10e3*1e-9\n",
    "print('Volume Change:',np.around(dVdt,1),'(km3/yr)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it! From just a few lines of code we can create a dataset that allows us to estimate the draw-down pattern using ICESat-2 for PIG for both space and time.   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Constructing a digital elevation model\n",
    "\n",
    "From **fittopo.py** we also have the mean surface height from all our solutions and which are now referenced to the same common time 2020. So, let’s go through a few steps to show how we can obtain a DEM in a relatively easy manner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_h = h_mod.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.scatter(lon[::100], lat[::100], s=3, c=h_h[::100], alpha=.7, transform=ccrs.PlateCarree(), cmap='terrain')\n",
    "plt.clim([0,1500])\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation (m)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let us remove the bad data that we identified in residual filtering!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h_h[np.isnan(z_i)] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(12,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.scatter(lon[::100], lat[::100], s=3, c=h_h[::100], alpha=.7, transform=ccrs.PlateCarree(), cmap='terrain')\n",
    "plt.clim([0,1500])\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation (m)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_h, y_h = transform_coord('4326', '3031', lon, lat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construct a grid for the DEM with a 1 km resolution using **make_grid** function or any other way you wish."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xdem, Ydem = make_grid(x_h.min(), x_h.max(), y_h.min(), y_h.max(), 1e3, 1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolate our grounded ice mask to our DEM-grid using bilinear interpoaltion (order=1). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Zm_i = interp2d(Xm, Ym, Zm, Xdem.ravel(), Ydem.ravel(), order=1).reshape(Xdem.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These interpolators are avaliable in **utils.py** and are simplfied versions of stand-alone programs that can be used for data analysis inside notebooks etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils import interpkrig, interpmed, interpgaus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Identify and remove any NaN's in the data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_nan = ~np.isnan(h_h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some notes on interpolation algorithms \n",
    "\n",
    "There are many different types of interpolation functions, which all have pros and cons, and some are more appropriate than others to specific problems. One thing to note is if you have very dense datasets with good point density and coverage there will usually be very little difference in the interpolated results from the different algorithms. One should note that to obtain good results data should be relatively clean and free of outliers!\n",
    "\n",
    "**Some common interpolation algorithm**\n",
    "- ***Kriging/Collocation (interpkrig):***\\\n",
    "Least-Squares Collocation (LSC), mostly known as Kriging, is a geostatistical method of spatial interpolation modeled by a gaussian process governed by the data covariance. In general the LSC method is better at preserving spatial detail and variability than other simpler distance-based interpolation methods, if the spatial auto-correlation of the data is substantial. LSC, also generates estimates of uncertainty for each prediction and also allows for error-propagation via the covariance structure.\n",
    "- ***Gaussian Weighed Average (interpgaus):***\\\n",
    "One of the most common interpolation approaches to interpolate irregular scattered data onto a regular grid using a weighted average, with a gaussian weighing function based on distance away from the prediction point. This interpolation approach provides a smooth continuous surface constrained by the range of the data. However, the weighted average is sensitive to large outliers and the data distribution, the same as for the median.\n",
    "- ***Median (interpmed):***\\\n",
    "Using the median for spatial interpolation has some advantages, as it is less sensitive to outliers in the data, and thus provide very clean and nice-looking grids. However, for producing volume change estimates for large regions (continental) it is less useful. This is because, as we could see in the histogram, we usually have a negative skewed distributions of elevation change. Thus, using the median we will estimates value will be biased in the positive direction. Further, the usage of the median does not allow us to incorporate far-field data, as no weighting for distance is possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z_dem = interpkrig(x_h[i_nan][::100], y_h[i_nan][::100], h_h[i_nan][::100].copy(),\\\n",
    "                   np.ones(h_h[i_nan][::100].shape)*1, Xdem, Ydem, n=20, d=20e3, a=5e3)[0]\n",
    "\"\"\"\n",
    "z_dem = interpgaus(x_h[i_nan][::100], y_h[i_nan][::100], h_h[i_nan][::100].copy(),\\\n",
    "                  np.ones(h_h[i_nan][::100].shape)*1, Xdem, Ydem, n=20, d=20e3, a=5e3)[0]\n",
    "                  \n",
    "z_dem = interpmed(x_h[i_nan][::100], y_h[i_nan][::100], h_h[i_nan][::100].copy(), Xdem, Ydem, n=20, d=20e3)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output of the functions are 1D-arrays, as you might want to interpolate points data as well (flight-tracks etc), so it needs to be changed back to a grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z_dem = z_dem.reshape(Xdem.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apply the mask to the grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Z_dem[Zm_i == 0] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(np.flipud(Z_dem), cmap='terrain', vmin=0, vmax=1400)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality assesment using crossover analysis\n",
    "\n",
    "One of the most widely used methods to study elevation change using altimetry data is the **\"crossover\"** method. It is the most robust and accurate method currently available. It can be used for both deriving elevation change and to study the precision of the measurements themselves.\n",
    "\n",
    "![crossover_method.png](crossover_method.png \"Figure of Eleavtion Change from ICESat\")\n",
    "\n",
    "\n",
    "From this method the surface elevation change is estimated at the intersection of two sub-satellite ground tracks at different times. At this crossover location the two closest points in each track is used to interpolate both elevation and time to the crossover location, allowing us to compute $\\delta h/\\delta t$ = $H_{t_{2}}-H_{t_{1}}$. [Image credit Zwally et al. (1989)]\n",
    "\n",
    "**\"captoolkit\"** has a program called **xover.py** to compute satellite/airborne crossover differences! Let us look at the input to the program!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python xover.py -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./xover.py ../../tutorial-data/gridding-time-series/org/IS2_HEIGHTS_PIG_A.h5 ../../tutorial-data/gridding-time-series/org/IS2_HEIGHTS_PIG_D.h5 \\\n",
    "-o ./xovers.h5 -r 100 -p 3031 -d 10 -k 20 20 -m linear -v orbit lon lat t_year h_elv dum dum dum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!h5ls ../../tutorial-data/gridding-time-series/proc/xovers.h5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('../../tutorial-data/gridding-time-series/proc/xovers.h5','r') as f_xo:\n",
    "    lon_xo = f_xo['lon'][:]\n",
    "    lat_xo = f_xo['lat'][:]\n",
    "    dh_xo  = f_xo['dh_elv'][:]\n",
    "    dt_xo  = f_xo['dt_year'][:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(9,10))\n",
    "ax = plt.axes(projection=ccrs.SouthPolarStereo())\n",
    "plt.scatter(lon_xo[np.abs(dt_xo)>6./12], lat_xo[np.abs(dt_xo)>6./12], s=5,\\\n",
    "            c=dh_xo[np.abs(dt_xo)>6./12]/dt_xo[np.abs(dt_xo)>6./12], \\\n",
    "            alpha=.7, transform=ccrs.PlateCarree(), cmap='coolwarm_r')\n",
    "plt.clim([-1.,1.])\n",
    "plt.colorbar(fraction=0.0320, pad=0.02, label='Elevation Change (m/yr)')\n",
    "ax.coastlines('50m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_xo, y_xo = transform_coord('4326','3031',lon_xo, lat_xo)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(dt_xo, dh_xo,'.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(dt_xo[np.abs(dh_xo)<5],dh_xo[np.abs(dh_xo)<5],'.',markersize=0.2)\n",
    "c_fit = np.polyfit(dt_xo[np.abs(dh_xo)<5], dh_xo[np.abs(dh_xo)<5], 1)\n",
    "plt.plot(dt_xo[np.abs(dh_xo)<5], np.polyval(c_fit, dt_xo[np.abs(dh_xo)<5]),'r')\n",
    "plt.ylabel('dh (m)')\n",
    "plt.xlabel('dt (yrs)')\n",
    "print('Trend:',np.around(c_fit[0],2),'(m/yr)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the quality assesment lets only use height differences of < 5 meters and time differences < 30 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i_keep = (np.abs(dh_xo) < 5) & (np.abs(dt_xo) < 1./12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_xo, y_xo, dt, dh = x_xo[i_keep], y_xo[i_keep], dt_xo[i_keep], dh_xo[i_keep] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.hist(dh[np.abs(dh) < 1.5], 100,edgecolor='k')\n",
    "plt.xlabel('dh less than 30 days (m)')\n",
    "plt.ylabel('Frequency')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('A/D-Bias:',np.around(np.nanmedian(dh)*100,2), 'cm')\n",
    "print('Error:',np.around(mad_std(dh)*100,2), 'cm')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's read some surface slope data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xs, Ys, Zs = tiffread('../../tutorial-data/gridding-time-series/REMA_1km_slope_filled_deg.tif')[0:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interpolate slope values from raster to crossovers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = interp2d(Xs, Ys, Zs, x_xo, y_xo, order=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bin data as a function of surface slope and output some statistics!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slp_b, bias_b, error_b = binning(slope, dh, xmin=0,xmax=1, dx=0.025, window=0.025)[0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(slp_b, bias_b, '-o')\n",
    "plt.ylabel('Bias (m)')\n",
    "plt.xlabel('Slope (deg)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(slp_b, error_b, '-o')\n",
    "plt.ylabel('Error (m)')\n",
    "plt.xlabel('Slope (deg)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope.min(), slope.max(), slope.mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Some homework maybe?\n",
    "\n",
    "If there is enough time let us breakout in different groups and try to work with the data in different ways:\n",
    "\n",
    "Try to play around with different parameters and see how they impact your results.\n",
    "- Search radius, grid resolution and model selection etc. in \"fittopo.py\".\n",
    "- Spatio-temporal resolution and correlation lengths in \"intergaus3d.py\".\n",
    "- Change the thresholds for the global and regional outlier editing"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
